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Abstract 

J*. 

Metamaterials are constructed such that, for a narrow range of frequencies, 
the momentum density depends on the local displacement gradient, and the stress 
depends on the local velocity. In these models the momentum density generally 
depends not only on the strain, but also on the local rotation, and the stress is 
generally not symmetric. A variant is constructed for which, at a fixed frequency, 
the momentum density is independent of the local rotation (but still depends on 
the strain) and the stress is symmetric (but still depends on the velocity). Gen- 
fSJ | eralizations of these metamaterials may be useful in the design of elastic cloaking 

devices. 
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q ; 1 Introduction 

The possibility of cloaking objects to make them invisible has attracted substantial recent 
attention. Alii and Engheta (p], [2]) following earlier work of Kerker [3], found that the 
scattering from an object could be substantially reduced by surrounding it by a plasmonic 
or metamaterial coating. Milton and Nicorovici [I], expanding on the work of [5] and [6], 
proved that the dipole moments of a polarizable point dipole or clusters of polarizable line 
dipoles would vanish when placed within a cloaking region surrounding a superlens, and 
this has been confirmed numerically [7]. Bruno And Lintner [8] show that only partial 
cloaking occurs when the polarizable object is not discrete. Ramm [9] and Miller pU] 
found that cloaking could be achieved by sensing and manipulating the fields near the 
boundary of the object. 

Perhaps the greatest interest has been generated by transformation based cloaking. 
This type of cloaking was first discovered by Greenleaf Lassas and Uhlmann ([11]. [T2] ) in 
the context of electrical conductivity. Their idea was to apply the well known fact (see, 
for example, [13]) that the equations of electrical conductivity retain their form under 
coordinate transformations, using a singular transformation which mapped a point to a 
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sphere (the cloaking region), and which was the identity outside a larger sphere. Perturb- 
ing the conductivity in the cloaking region corresponds to changing the conductivity at 
a point in the equivalent problem, and this has no effect on the fields outside the larger 
sphere. This type of cloaking was extended to the full time harmonic Maxwell equa- 
tions by Pendry, Schurig and Smith [T3] , and at the same time Leonhardt [15] discovered 
a different transformation based, two-dimensional, geometric optics or high frequency 
acoustic cloaking scheme which only utilized isotropic materials. The former cloaking 
has been verified numerically ([IE]; p2]; [IS]), placed on a firm theoretical foundation 
([19]), and experimentally substantiated in the microwave regime in an approximate way 
in two- dimensions by Schurig, Mock, Justice, Cummer, Pendry, Starr, and Smith [20] us- 
ing Schelkunoff and Friis's idea (plj) of utilizing metamaterials with split ring resonators 
to achieve artificial magnetism. Cai, Chettiar, Kildishev, and Shalaev [22] proposed an- 
other design which may experimentally achieve approximate two-dimensional cloaking in 
the visible regime. 

Milton, Briane and Willis [23] found that applying transformation based cloaking to 
continuum elastodynamics would require new materials with very unusual properties. 
This is because the continuum elastodynamic equations do not retain their form under 
coordinate transformations, but rather take the form of the equations that Willis ([24], 
[25], [26]) found described the ensemble averaged behavior of composite materials. One 
interesting feature is that the density at a given frequency is required to be matrix valued. 
Willis [27] had proved that his density operator had this property. It also follows from 
the work of Movchan and Guenneau [28J and Avila, Griso, and Miara [29], that there 
exist high contrast materials with a local anisotropic effective density. Simple models 
with anisotropic and possibly complex density were obtained by Milton and Willis [30J. 
These microstructures generalized a construction of Sheng, Zhang, Liu, and Chan [31] 
and Liu, Chan, and Sheng [32j, who established that the density can be negative over a 
range of frequencies. The same conclusion was reached and rigorously proved by Avila, 
Griso, and Miara [25]. A comparison of Maxwell's equations, which can be written in 
the form 

° r C m ^)={J l eE} J , (1.1) 



dxi \ dxk 
where 

[in which E is the electric field, e the electric permittivity tensor, \x the magnetic perme- 
ability tensor, and e ijm = 1 (-1) if ijm is an even (odd) permutation of 123 and is zero 
otherwise] with the equations of continuum elastodynamics 

TT~ ^ijki— = -{to pujj, (1.3) 



O I % IKK, r\ 

OXi \ OXk 

[in which u is the displacement field, p is the density, and C is now the elasticity tensor] 
suggests that p as a function of the frequency to could share many of the same proper- 
ties as e(uj). This has been established, and in fact for every function p(uj) satisfying 
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these properties one can construct a model which has approximately that function as its 
effective density as a function of frequency ([3D]). Cummer and Schurig [33J investigated 
transformation based acoustic cloaking in two- dimensions and, due to a mathematical 
equivalence with the electromagnetic problem, found that it could be achieved provided 
one could construct the necessary materials with anisotropic density. 

Anisotropic density is not the only new property needed to achieve transformation 
based elasticity cloaking. One needs materials where the constitutive law, like that in 
the Willis equations, couples the stress with not only with the strain but also with the 
velocity, and couples the momentum density not only with the velocity but also with the 
displacement gradient through the strain. Also, unlike in the Willis equations, one wants 
this dependence to be local in space and to apply to a single microgeometry rather than 
to an ensemble average of microgeometries. Here we show that such unusual behavior 
can be realized in a model such that, for a narrow range of frequencies, the associated 
waves have wavelength much larger than the microstructure, and the momentum density 
depends on the displacement gradient, and the stress depends on the velocity. Curiously 
the momentum density depends not only on the strain, but also on the local rotation, 
and the stress is not symmetric. We also construct a variant of the model for which (at 
a fixed frequency) the momentum density is independent of the local rotation (but still 
depends on the strain) and the stress is symmetric (but still depends on the velocity). 
Both models rely heavily on the discovery Sheng, Zhang, Liu, and Chan [31] that one 
can design structures which, at a fixed frequency, respond as if they have negative mass. 
Following their ideas, and the simplified constructions of Milton and Willis [30] which 
incorporate these ideas, figure 1 illustrates a structure with negative mass. In the models 
considered in this paper both positive and negative masses are idealized as point masses. 

This work on metamaterials with macroscopic behavior outside that of continuum 
elastodynamics, even though they are governed by continuum elastodynamics at the mi- 
croscale, is preceeded by work on metamaterials with macroscopic non-Ohmic, possibly 
non-local, conducting behavior, even though they conform to Ohm's law at the mi- 
croscale, ([31]; [35]; [36]; [37]; [38]; [39]) by work on metamaterials with non-Maxwellian 
macroscopic electromagnetic behavior, even though they conform to Maxwell's equations 
at the microscale ([40]), and by work on metamaterials with a macroscopic higher or- 
der gradient or non-local elastic response even though they are governed by usual linear 
elasticity equations at the microscale ([S]; [12]; [43]). 

For simplicity we consider a two dimensional model, although it can be generalized 
to three dimensions. The model is illustrated in figure 2. 

2 The momentum density in the model 

Figure 3 shows the unit cell in the model. The masses are approximated as point masses, 
and the springs may be taken to have all the same spring constant hK, which scales in 
proportion to the size of the unit cell. This ensures that the spring network, by itself, 
responds elastic material in the limit h — > 0. 



3 



Figure 1: This structure responds as if it had negative mass to oscillations at a fixed 
frequency above resonance, such that the spherical core oscillates out of phase with the 
rigid, but light, surrounding shell. The springs have the same spring constant to ensure 
that the effective mass is isotropic. 




Figure 2: The model with strange elastic behavior. The large solid black circles have 
positive mass, while the neighbouring large white circles have negative mass at the given 
frequency. They are connected to the spring network by rods of fixed length, which 
alternatively can be regarded as springs with infinite spring constants. 
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Figure 3: The unit cell of the model. 



Without loss of generality let us place the origin at the center x of the unit cell under 
consideration. Then when the material is a rest the points A, B, C, D, E and F in figure 
3 are at positions 



x A = (-h,0), x B = (0, -h), 
X£ = (— ch, 0), Xi? = (c/i, 0), 



x c = (M), x D = (0,/i) 



(2-1) 



where c is a parameter between and 1 which controls the inclination of the rods joining 
E and F to B and D: these rods have length hy/1 + c 2 . In the limit ft^Owe assume 
the (infinitesimal) physical displacements at the points x^, x#, and X£> derive from 
some smooth complex valued displacement field u(x), and in particular 



u A = Ke[u(x A )e-^] 
u B = 3te[u(x B )e- iwt ] 
u c = 3fi e [u(x A )e- <wt ] 
u D = Ke[u(x D )e-^] 



sfte[(u -/iq)e-^], 
3?e[(u - W)e-i 
3?e[(u + /iq)e-^], 
Ke[(u + H^'l 



where 



u = u(x ), q 



<9u 

dxi 



w 



X=XQ 



<9u 

dx 2 



(2.2) 



(2.3) 



X=XQ 
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and xo = because the origin is at the center of the unit cell. 

At the macroscopic level we choose not keep track of the displacements u E and u F . 
These are internal hidden variables which however can be recovered from u# and Up. 
They have the form 

u E « 3?e[(u + hs)e- iuit }, 

u F « Ke[(u - hs)e' iult }, (2.4) 

where the complex vector s remains to be determined. Since there is a rigid rod connect- 
ing the points D and E we have, in the limit h — *■ 0, the constraint that w — s must be 
perpendicular to the line joining D and E, i.e. 

= (w - s) • (c, 1) = c(wx - s x ) + (w 2 - s 2 ). (2.5) 

Similarly since there is a rigid rod connecting the points B and E we have, in the limit 
h — > 0, the constraint 

= (w + s) ■ (-c, 1) = -c(wi + si) + (w 2 + s 2 ). (2.6) 

These equations have the solution 

Sl = ty 2 /c = — — , s 2 =CWi = C-— . (2.7) 

c ax 2 ax 2 

One can easily see that if c is replaced by — c then s gets replaced by — s, which justifies 
the formula ( 12. 4p for up. 

Now suppose the masses at the points E and F, which we call a mass pair, are 
respectively chosen to have the form 

rriE = hm, m E = —hm + Sh 2 , (2-8) 

where m is a positive real constant, and 8 is a possibly complex constant with a non- 
negative imaginary part. Then at any time t the physical momentum in the unit cell is, 
to second order in h, 

ne{-iLo[m E (uo + hs) + m F (u - hs)]e- iuJt } w h 2 ^e{-iuj[2ms + 5u ]e' iuJt }. (2.9) 

Since the area of the unit cell is 2h 2 we see that the complex momentum density is 

p = — icums + (<5/2)(— iuuo), (2-10) 

in which — iuuQ is the complex velocity of the unit cell, and s depends on the deformation 
gradient through (12.71) . 
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3 The stress in the model 



The acceleration of the material will also generate stress. To see this, note that to leading 
order in h, the masses at E and F must be accelerated by complex forces F = —u 2 mhu 
and — F acting on the respective masses (the physical forces are obtained by multiplying 
these by e~ luJt and taking the real part). Let Fed and Fj?d be the forces which the rods 
ED and FD exert on the vertex D and let Feb and Feb be the forces which the rods 
EB and FB exert on the vertex B. Balance of forces at E and F requires that 

Fed + Feb = — F, Fed + Feb = F. (3-1) 

Since these forces are aligned with their respective rods we have 

F E D = -F F B = a(c,l), F EB = -F F D=P(c,-l), (3.2) 

where the constants a and /3 are such that ( 13 .ip is satisfied: 

(a + @)c = -Fx, a- (3 = -F 2 . (3.3) 

Forgetting for a moment the springs and surrounding material, to maintain the motions 
a force 

ht = -Fed ~ Fed = -a(c, 1) + P(c, -1) = (cF 2 , Fjc) = -cu 2 mh(cu 2 , Ul /c) (3.4) 

needs to act on the vertex D to balance the forces which the rods ED and FD exert 
on this vertex. Similarly a force —ht needs to act on the vertex B to balance the forces 
which the rods EB and FB exert on this vertex. On the other hand, no additional forces 
need to be exerted on the vertices A and C . 

Now consider a rectangular sample (with dimensions, say of order y/h, small com- 
pared to the wavelength) as in figure 2 with (complex) forces acting on the boundary 
vertices to maintain the motions. These forces will have two components: an elastic (or 
possibly viscoelastic) component to compensate for the (complex) stress cte caused by 
the (complex) strain in the spring network, and an inertial component to compensate for 
the inertial stress crj caused by the acceleration — u 2 u of the material. The compensating 
inertial component will be zero on the left and right sides of the network. At the top 
5 vertices in figure 2 the compensating inertial component will be approximately 2ht at 
each vertex. The extra factor of 2 arises because of the extra load that the springs at 
the top boundary carry to support the inertial component of the forces at the 4 vertices 
immediately below the top. (At the other interior nodes these forces are balanced). Thus 
the compensating inertial component per unit length, is t on the top, — t on the bottom, 
and zero on the sides. By the definition of stress this should be equated with cr^n, where 
n is the unit normal to the boundary. Thus we deduce that 

<"=(n t) = (n ~i mCU j). (3.5) 
\Q h) \0 -id mui/c J y ' 



7 



Of course, in the limit h — > 0, the stress in the spring network (excluding the masses and 
connecting rods) will be governed by a normal relation 

(T E = C[Vu+(Vuf], (3.6) 

where C is the elasticity tensor of the network. 

Our choice to define the macroscopic stress through the forces at the vertices of the 
unit cell (rather than as a volume average of the microscopic stress field, which would 
result in a symmetric stress) makes good physical sense and is consistent with our choice 
to define the macroscopic displacement field through the values of the displacement at 
the vertices of the unit cell. These choices bear some resemblance to the way Pendry, 
Holden, Robbins and Willis [Hj define macroscopic electromagnetic fields through field 
values at the boundary of the unit cell. In defining the stress in this way it is important 
that the unit cell be chosen so that each mass pair (of positive and negative masses) is 
not split by the boundary of the unit cell: each mass pair and their four connecting rods 
is regarded as an indivisible unit in the same way that in defining the electric polarization 
field one would not choose to place the boundary in a dielectric material so as to split 
the positive and negative charges of the constituent atoms. 



4 Summary 



In the limit h — > the constitutive law takes the form 



C S 
D p 



Vu 
v 



(4.1) 



where a = cte + &i is the total stress, v = — ium is the velocity, and p = (5/2)1 is the 
density. In the basis where cr, p, Vu and v are represented by the vectors 



C21 

\CT22 / 



P 



Pi 
P2 



Vu 



( dui/dxi \ 
du2/dxi 
dui/dx 2 

\du 2 /dx 2 ) 



Vl 

v 2 



(4.2) 



the third-order tensors S and D, from (12. 7p . (12.101) and (I3.5p . are represented by the 
matrices 



/ 









V —iujm/c 






-iujmc 




D 



—ium/c 
—iumc 



(4.3) 



Thus the matrix entering the constitutive law (14.11) is symmetric. Since the constitutive 
law and the equation of motion are asymptotically independent of h, the waves associated 
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2h 




Figure 4: The unit cell of a generalized model. 

with their solutions will have wavelength much larger that the microstructure. In this 
circumstance the usual continuum elastodynamic equations normally apply. However 
this model shows that this is not always the case. 

The model has some unusual features. Although the net amount of mass in a region 
of unit area is independent of h, the total amount of positive (or negative) mass in a 
region of unit area scales like 1/h. However in any physical model, as usual, one never in 
fact takes the limit h — > 0, but instead one sets the microstructure as small as required 
to get a reasonable approximation to the limiting behavior for a desired set of applied 
fields. Also the amount of positive (or negative) mass in the unit cell only represents 
the amount of apparent mass, not gravitational mass, which could be quite different if 
the apparent mass arises from substructures which are close to resonance. Of course 
the analysis given here needs some rigorous justification to check, for example, that the 
displacements can be approximated by a smooth field u(x) in the limit h — > 0. Also from 
a practical point of view it needs to be checked that the behavior is stable to sufficiently 
small variations in the masses or spring constants from cell to cell. 

The model can easily be generalized. One simple generalization has additional masses 
and connecting rods as illustrated in figure 4, with masses 

mo = —hm + 5h 2 , itlh = hm, (4.4) 

at the points G and H which, when the material is at rest, are located at 

x G = (-c7i,0), x H = (c7i,0), (4.5) 
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where d and ml are positive constants. By superposition the constitutive law will take the 
form (14.11) with a density p = 51 and with the third-order tensors S and D represented 
by the matrices 



\ 








iuim'd 



mc) 



D 



\ioj(m!/d — m/c) 



iujim'd 







/ 



mc) 



iu){m! jd — m/c) 




(4.6) 



In particular if, at a given frequency, m'd = mc we obtain a model in which the mo- 
mentum does not depend on the local rotation and the stress is symmetric. With this 
condition satisfied one has the relations 



<r = CVu + S'u, V • cr = -iup = D'Vu - to 2 pu, (4.7) 

where 

S' = -iuS, D' = -zcuD, (4.8) 
and the only non-zero Cartesian elements of these two tensors are from (j4.ip . (14.21) and 

^22i = D' 122 = uj\m!jd - m/c). (4.9) 

The general form of these tensors S' and D' corresponds with that required for elasticity 
cloaking: they are real and satisfy the identities S'^ k = Sj ik = D' ki -. They may serve 
as a good starting point in the quest to find materials with a suitable combination of 
properties C, S', D' and p needed for elasticity cloaking. However, quite apart from 
this, the novel properties of this new class of metamaterials may have other important 
applications. 

In the design of these models perhaps the key feature is that in the narrow frequency 
range under consideration the amount of negative mass in the unit cell (due to hidden 
internal masses moving out of phase with the motion) almost balances the amount of 
positive mass in the unit cell. This ensures that the momentum density approaches a 
constant as the cell size approaches zero, while keeping non-trivial the stresses caused by 
the moving masses in the unit cell. It will be interesting to see if other models with this 
key feature also have a local constitutive law of the form (14.11) . with non-zero coupling 
terms S and D. 
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